********************************************************************************
* Competition model with bi-level model and disjunction, using convex hull and 
* KKT conditions to solve the model.
* Lei Zhang
* 10/22/2013
********************************************************************************

SETS
    i   our_facilities          /F1/
    ii  competitor's_facilities /FC1/
    j   markets                 /M1*M2/
    k   products                /LOX/
    t   time_periods            /T1*T2/;

SCALAR
    CMAX    upper_bound_of_capacity_expansion   /100.0/
    U       upper_bound_of_capacity             /500.0/
    DCC     upper_bound_of_competitor_expansion /100.0/
    M       Big-M_method                        /1000.0/
    delta   Coefficient_of_competitor_expansion /1.1/;

TABLE falpha(i,k,t) Fixed_cost_to_expand_capacity
            T1      T2
F1.LOX      10.0    9.0;

TABLE fbeta(i,k,t)  Unit_investment_cost_to_expand_capacity
            T1      T2
F1.LOX      1.1     0.9;

TABLE fgamma(i,k,t) Unit_production_cost
            T1      T2
F1.LOX      0.6     0.6;

TABLE TR(i,j,k) Unit_transportation_cost
            LOX
F1.M1       0.1
F1.M2       0.2;

TABLE D(j,k,t)  Demand_of_market
            T1      T2
M1.LOX      35.0    38.0
M2.LOX      20.0    25.0;

TABLE P(i,j,k,t)    Price
            T1      T2
F1.M1.LOX   2.0     2.0
F1.M2.LOX   2.0     2.0;

TABLE PC(ii,j,k,t)  Competitor's_price
            T1      T2
FC1.M1.LOX  2.1     2.1
FC1.M2.LOX  2.1     2.1;

*------------------------------------------------------------------------------*

VARIABLES
    npv             Net_present_value
    lambda1(j,k,t)  KKT;

BINARY VARIABLES
    x(i,k,t)        Selection_of_capacity_expansion

    w1(i,k,t)       Big-M_KKT
    w2(ii,k,t)      Big-M_KKT
    w3(k,t)         Big-M_KKT
    w4(k,t)         Big-M_KKT
    w5(k,t)         Big-M_KKT
    w6(ii,k,t)      Big-M_KKT
    w7(ii,k,t)      Big-M_KKT
    w8(i,k,t)       Big-M_KKT
    w9(i,k,t)       Big-M_KKT
    w10(ii,k,t)     Big-M_KKT
    w11(ii,k,t)     Big-M_KKT;

POSITIVE VARIABLES
    y(i,j,k,t)      Amount_of_product_that_facility_sells_to_market
    yc(ii,j,k,t)    Amount_of_product_that_competitor's_facility_sells_to_market
    c(i,k,t)        Capacity_of_facility
    dc(i,k,t)       Capacity_expansion
    cc(ii,k,t)      Capacity_of_competitor's_facility
    c1(i,k,t)       Convex_hull
    c2(i,k,t)       Convex_hull
    cc1(ii,k,t)     Convex_hull
    cc2(ii,k,t)     Convex_hull
    z(k,t)          Convex_hull

    mu1(i,k,t)      KKT
    mu2(ii,k,t)     KKT
    mu3(k,t)        KKT
    mu4(k,t)        KKT
    mu5(k,t)        KKT
    mu6(ii,k,t)     KKT
    mu7(ii,k,t)     KKT
    mu8(i,k,t)      KKT
    mu9(i,k,t)      KKT
    mu10(ii,k,t)    KKT
    mu11(ii,k,t)    KKT;

*------------------------------------------------------------------------------*

EQUATIONS
    obj             Objective_function_maximize_the_NPV
    epd(i,k,t)      Investment_decision_in_capacity_expansion
    upb(i,k,t)      Upper_bound_of_capacity_expansion
    stn1(i,j,k,t)   KKT_stationary
    stn2(ii,j,k,t)  KKT_stationary
    stn3(k,t)       KKT_stationary
    tmd(j,k,t)      Demand_satisfaction_for_all_markets
    cvx1(i,k,t)     Convex_hull
    cvx2(ii,k,t)    Convex_hull
    cpd1(ii,k,t)    Competitor_expansion_disjunction_1
    dmd(i,k,t)      Demand_satisfaction_to_the_installed_capacity
    dmdc(ii,k,t)    Demand_satisfaction_to_the_competitor
    cvx(k,t)        Convex_hull_variable
    cvx3(k,t)       Convex_hull
    cvx4(k,t)       Convex_hull
    cvx5(ii,k,t)    Convex_hull
    cvx6(ii,k,t)    Convex_hull
    cvu1(i,k,t)     Convex_hull
    cvu2(i,k,t)     Convex_hull
    cvu3(ii,k,t)    Convex_hull
    cvu4(ii,k,t)    Convex_hull
    cpt11(i,k,t)    Big-M_KKT_complementary_slackness
    cpt12(i,k,t)    Big-M_KKT_complementary_slackness
    cpt21(ii,k,t)   Big-M_KKT_complementary_slackness
    cpt22(ii,k,t)   Big-M_KKT_complementary_slackness
    cpt31(k,t)      Big-M_KKT_complementary_slackness
    cpt32(k,t)      Big-M_KKT_complementary_slackness
    cpt41(k,t)      Big-M_KKT_complementary_slackness
    cpt42(k,t)      Big-M_KKT_complementary_slackness
    cpt51(k,t)      Big-M_KKT_complementary_slackness
    cpt52(k,t)      Big-M_KKT_complementary_slackness
    cpt61(ii,k,t)   Big-M_KKT_complementary_slackness
    cpt62(ii,k,t)   Big-M_KKT_complementary_slackness
    cpt71(ii,k,t)   Big-M_KKT_complementary_slackness
    cpt72(ii,k,t)   Big-M_KKT_complementary_slackness
    cpt81(i,k,t)    Big-M_KKT_complementary_slackness
    cpt82(i,k,t)    Big-M_KKT_complementary_slackness
    cpt91(i,k,t)    Big-M_KKT_complementary_slackness
    cpt92(i,k,t)    Big-M_KKT_complementary_slackness
    cpt101(ii,k,t)  Big-M_KKT_complementary_slackness
    cpt102(ii,k,t)  Big-M_KKT_complementary_slackness
    cpt111(ii,k,t)  Big-M_KKT_complementary_slackness
    cpt112(ii,k,t)  Big-M_KKT_complementary_slackness;
        
obj..               npv =e= SUM((i,j,k,t), P(i,j,k,t)*y(i,j,k,t))
                        - SUM((i,k,t), falpha(i,k,t)*x(i,k,t))
                        - SUM((i,k,t), fbeta(i,k,t)*dc(i,k,t))
                        - SUM((i,k,t), fgamma(i,k,t)*c(i,k,t))
                        - SUM((i,j,k,t), TR(i,j,k)*y(i,j,k,t));
epd(i,k,t)..        c(i,k,t) =e= c(i,k,t-1) + dc(i,k,t);
upb(i,k,t)..        dc(i,k,t) =l= CMAX * x(i,k,t);
stn1(i,j,k,t)..     P(i,j,k,t) + lambda1(j,k,t) + mu1(i,k,t) =e= 0.0;
stn2(ii,j,k,t)..    PC(ii,j,k,t) + lambda1(j,k,t) + mu2(ii,k,t) =e= 0.0;
stn3(k,t)..         mu3(k,t) + delta*SUM(j, D(j,k,t))*mu4(k,t)
                        + delta*SUM(j, D(j,k,t))*mu5(k,t)
                        + SUM(ii, DCC*mu6(ii,k,t))
                        + SUM(i, U*mu9(i,k,t)-U*mu8(i,k,t))
                        + SUM(ii, U*mu11(ii,k,t)-U*mu10(ii,k,t)) =e= 0.0;
tmd(j,k,t)..        SUM(i, y(i,j,k,t)) + SUM(ii, yc(ii,j,k,t)) 
                        - D(j,k,t) =e= 0.0;
cvx1(i,k,t)..       c(i,k,t) =e= c1(i,k,t) + c2(i,k,t);
cvx2(ii,k,t)..      cc(ii,k,t) =e= cc1(ii,k,t) + cc2(ii,k,t);
cpd1(ii,k,t)..      cc1(ii,k,t) =e= cc1(ii,k,t-1);
dmd(i,k,t)..        SUM(j, y(i,j,k,t)) - c(i,k,t) =l= 0.0;
dmdc(ii,k,t)..      SUM(j, yc(ii,j,k,t)) - cc(ii,k,t) =l= 0.0;
cvx(k,t)..          z(k,t) - 1.0 =l= 0.0;
cvx3(k,t)..         delta*SUM(j, D(j,k,t))*z(k,t) - SUM(i, c1(i,k,t))
                        - SUM(ii, cc1(ii,k,t)) =l= 0.0;
cvx4(k,t)..         SUM(i, c2(i,k,t)) + SUM(ii, cc2(ii,k,t))
                        - delta*SUM(j, D(j,k,t))*(1.0-z(k,t)) =l= 0.0;
cvx5(ii,k,t)..      cc2(ii,k,t) =l= cc2(ii,k,t-1) + DCC*(1.0-z(k,t));
cvx6(ii,k,t)..      cc2(ii,k,t-1) - cc2(ii,k,t) =l= 0.0;
cvu1(i,k,t)..       c1(i,k,t) =l= U * z(k,t);
cvu2(i,k,t)..       c2(i,k,t) =l= U * z(k,t);
cvu3(ii,k,t)..      cc1(ii,k,t) =l= U * (1.0-z(k,t));
cvu4(ii,k,t)..      cc2(ii,k,t) =l= U * (1.0-z(k,t));
cpt11(i,k,t)..      mu1(i,k,t) =l= M * (1.0-w1(i,k,t));
cpt12(i,k,t)..      SUM(j, y(i,j,k,t)) - c(i,k,t) =g= -M * w1(i,k,t);
cpt21(ii,k,t)..     mu2(ii,k,t) =l= M * (1.0-w2(ii,k,t));
cpt22(ii,k,t)..     SUM(j, yc(ii,j,k,t)) - cc(ii,k,t) =g= -M * w2(ii,k,t);
cpt31(k,t)..        mu3(k,t) =l= M * (1.0-w3(k,t));
cpt32(k,t)..        z(k,t) - 1.0 =g= -M * w3(k,t);
cpt41(k,t)..        mu4(k,t) =l= M * (1.0-w4(k,t));
cpt42(k,t)..        delta*SUM(j, D(j,k,t))*z(k,t) - SUM(i, c1(i,k,t))
                        - SUM(ii, cc1(ii,k,t)) =g= -M * w4(k,t);
cpt51(k,t)..        mu5(k,t) =l= M * (1.0-w5(k,t));
cpt52(k,t)..        SUM(i, c2(i,k,t)) + SUM(ii, cc2(ii,k,t))
                        - delta*SUM(j, D(j,k,t))*(1-z(k,t)) =g= -M * w5(k,t);
cpt61(ii,k,t)..     mu6(ii,k,t) =l= M * (1.0-w6(ii,k,t));
cpt62(ii,k,t)..     cc2(ii,k,t) - cc2(ii,k,t-1) - DCC*(1.0-z(k,t)) 
                        =g= -M * w6(ii,k,t);
cpt71(ii,k,t)..     mu7(ii,k,t) =l= M * (1.0-w7(ii,k,t));
cpt72(ii,k,t)..     cc2(ii,k,t-1) - cc2(ii,k,t) =g= -M * w7(ii,k,t);
cpt81(i,k,t)..      mu8(i,k,t) =l= M * (1.0-w8(i,k,t));
cpt82(i,k,t)..      c1(i,k,t) - U*z(k,t) =g= -M * w8(i,k,t);
cpt91(i,k,t)..      mu9(i,k,t) =l= M * (1.0-w9(i,k,t));
cpt92(i,k,t)..      c2(i,k,t) - U*(1.0-z(k,t)) =g= -M * w9(i,k,t);
cpt101(ii,k,t)..    mu10(ii,k,t) =l= M * (1.0-w10(ii,k,t));
cpt102(ii,k,t)..    cc1(ii,k,t) - U*z(k,t) =g= -M * w10(ii,k,t);
cpt111(ii,k,t)..    mu11(ii,k,t) =l= M * (1.0-w11(ii,k,t));
cpt112(ii,k,t)..    cc2(ii,k,t) - U*(1.0-z(k,t)) =g= -M * w11(ii,k,t);

*------------------------------------------------------------------------------*

MODEL COMPETITION /ALL/;
SOLVE COMPETITION USING MIP MAXIMIZING npv;

********************************************************************************
